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Abstract 

The 3-D nonequilibrium seeded air flow in the NAS A- 
Ames experimental MHD channel has been numeri- 
cally simulated. The channel contains a nozzle sec- 
tion, a center section, and an accelerator section 
where magnetic and electric fields can be imposed 
on the flow. In recent tests, velocity increases of up 
to 40% have been achieved in the accelerator section. 
The flow in the channel is numerically computed us- 
ing a 3-D parabolized Navier-Stokes (PNS) algorithm 
that has been developed to efficiently compute MHD 
flows in the low magnetic Reynolds number regime. 
The MHD effects are modeled by introducing source 
terms into the PNS equations which can then be 
solved in a very efficient manner. The algorithm has 
been extended in the present study to account for 
nonequilibrium seeded air flows. The electrical con- 
ductivity of the flow is determined using the program 
of Park. The new algorithm has been used to com- 
pute two test cases that match the experimental con- 
ditions. In both cases, magnetic and electric fields 
are applied to the seeded flow. The computed results 
are in good agreement with the experimental data. 

Introduction 

Magnetohyrodynamics (MHD) can be utilized to im- 
prove performance and extend the operational range 
of many systems. Potential applications include hy- 
personic cruise, advanced Earth-to-orbit propulsion, 
chemical and nuclear space propulsion, regenerative 
aerobraking, onboard flow control systems, test fa- 

* Graduate Research Assistant, Student Member AIAA 

Manager, Computational Fluid Dynamics Center, and Pro- 
fessor, Dept, of Aerospace Engineering. Fellow AIAA 

* Division Scientist, Associate Fellow, AIAA 
Copyright ©2005 by the American Institute of Aeronau- 
tics and Astronautics, Inc., all rights reserved. 


cilities, launch assist, and power generation. One of 
the critical technologies associated with these appli- 
cations is MHD acceleration. In order to study MHD 
acceleration, an experimental MHD channel has been 
built at NASA Ames Research Center by D. W. Bog- 
danoff, C. Park, and U. B. Mehta [1,2]. The channel 
is about a half meter long and contains a nozzle sec- 
tion, a center section, and an accelerator section. The 
channel has a uniform width of 2.03 cm. Magnetic 
and electric fields can be imposed upon the flow in 
the accelerator section. A cross section of the MHD 
channel is shown in Fig. 1. 

In the present study, the flow in the experimen- 
tal MHD channel is numerically simulated. Flow- 
fields involving MHD effects have typically been 
computed [3-15] by solving the complete Navier- 
Stokes (N-S) equations for fluid flow in conjunction 
with Maxwell’s equations of electromagnetodynam- 
ics. When chemistry and turbulence effects are also 
included, the computational effort required to solve 
the resulting coupled system of partial differential 
equations is extremely formidable. One possible rem- 
edy to this problem is to use the parabolized Navier- 
Stokes (PNS) equations in place of the N-S equa- 
tions. The PNS equations can be used to compute 
three-dimensional, supersonic viscous flowfields in a 
very efficient manner [16]. This efficiency is achieved 
because the equations can be solved using a space- 
marching technique as opposed to the time-marching 
technique that is normally employed for the complete 
N-S equations. 

Recently, the present authors have developed PNS 
codes to solve 2-D and 3-D supersonic MHD flowfields 
in both the high and low magnetic Reynolds number 
regimes [17-19]. The magnetic Reynolds number is 
defined as Re m = a e ^VooL where a e is the electrical 
conductivity, /x e is the magnetic permeability, is 
the ffeestream velocity, and L is the reference length. 
The new MHD PNS codes are based on NASA’s up- 
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wind PNS (UPS) code which was originally developed 
by Lawrence et al. [ 20 ]. The UPS code solves the PNS 
equations using Roe’s scheme in a fully conservative, 
finite- volume approach in general nonorthogonal co- 
ordinates. 

For many aerospace applications, including the 
present experimental MHD channel, the electrical 
conductivity of the fluid is low and hence the mag- 
netic Reynolds number is small. In these cases, it 
makes sense to use the low magnetic Reynolds num- 
ber assumption and reduce the complexity of the gov- 
erning equations. The MHD effects are modeled with 
the introduction of source terms into the PNS equar 
tions. 

Previously [19], the present authors used the low 
magnetic Reynolds PNS code to compute both 2 - 
D and 3-D flows in the NASA- Ames MHD channe l. 
These perfect gas (7 = 1.25) calculations assumed 
that the magnetic and electric fields, as well as the 
electrical conductivity, were constant in the accelera- 
tor section. In the present study, the 3-D simulations 
have been extended to include both equilibrium air 
flows as well as nonequilibrium seeded air flows. For 
the latter case, the electrical conductivity is variable 
and is computed using the program of Park [21]. 


Governing Equations 

Magnetogasdynamic Equations 

The governing equations for a viscous MHD flow 
with a small magnetic Reynolds number are given 
by [14]: 

Continuity equation 

§+V-(pV )=0 ( 1 ) 

Momentum equation 

^ + V- [pW + plj =V-f.+ JxB ( 2 ) 


The governing magnetogasdynamic equations are 
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where the superscript * refers to the nondimen- 
sional quantities. For convenience, the asterisks are 
dropped in the following equations. 

The governing equations written in vector form in 
a 3-D Cartesian coordinate system become 


dU dE* ( dF i t dG i dE v dF v dG v B 
~dT + = ~dx + ~dy’ + ~dT +SuHD 

( 6 ) 

where U is the vector of dependent variables, E it F* 
and G» are the inviscid flux vectors, and E v , F v and 
G v are the viscous flux vectors. The source term 
Smhd contains all of the MHD effects. The flux vec- 
tors are given by 
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Energy equation 


0 


^~-+^-[(pe t +p) V] = V-(V-f)— V-U+E-J (3) 
Ohm’s law 

J = <7 e (E + VxB) (4) 

where V is the velocity vector, B is the magnetic 
field vector, E is the electric field vector, and J is the 
conduction current density. 
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where 
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SmHD = Re-m 
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and the nondimensional shear stresses and heat fluxes 
are defined in the usual manner [16]. 

The governing equations are transformed into com- 
putational space and written in a generalized coordi- 
nate system (f,7j,C) as 


jU t + E e + F„ + G c = ^p ( 14 ) 

where 

e = (ky Ei . Ev)+ (ky Ti . Fv) 
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and J is the Jacobian of the transformation. 

The governing equations are parabolized by drop- 
ping the time derivative term and the streamwise 
direction (f) viscous flow terms in the flux vectors. 
Equation (14) can then be rewritten as 


Ef + 4 


Smhd 
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(16) 
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The primes in the preceding equations indicate that 
the streamwise viscous flow terms have been dropped. 

For turbulent flows, the two-layer B aldwin- Lomax 
turbulence model [22] has been modified to account 
for MHD effects. Only the expression for turbulent 
viscosity in the inner layer is changed. This modifi- 
cation for MHD flows is due to Lykoudis [23]. 

In order to “close” the preceding system of PNS 
equations, relations between the thermodynamic 
variables are required along with expressions for the 
transport properties /x and k. For a perfect gas, the 
pressure is computed from the relation 

p = (j-l)pe (18) 

where 7 = 7oo, and the transport properties are com- 
puted using Sutherland’s formulas [16]. For equilib- 
rium air computations, 7 and all other thermody- 
namic and transport properties are obtained from the 
simplified curve fits of Srinivasan et al. [24,25]. For 
nonequilibrium computations, the thermodynamic 
and transport properties are determined using the 
procedures described in the next section. 


Nonequilibrium Flow Equations 


For nonequilibrium flows, the species continuity equa- 
tions must be solved in addition to the magnetogas- 
dynamic equations given previously The magnet o- 
gasdynamic equations remain the same except for the 
additional term in the energy equation, which is due 
to the diffusion of the species. The nondimensional 
species continuity equations, expressed in 2-D trans- 
formed coordinates for a steady flow, are given by 


’ dc s f v \ dc s 1 . t 
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Re 00 dr ] 
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where c s is the mass fraction of species s, u s is the 
nondimensional production term, D is the nondimen- 
sional binary diffusion coefficient, and 0z = P ^ D °° . 
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The chemical model used in the present calcula- 
tions is similar to the clean-air model of Blottner et 
al. [26] and Prabhu et al. [27]. It consists of molecular 
oxygen (0 2 ), atomic oxygen (O), molecular nitrogen 
(iV 2 ), atomic nitrogen (N), nitric oxide {NO), ni- 
tric oxide ion (NO+) and electrons (e“). The follow- 
ing reactions are considered between the constituent 
species. 

(1) 0 2 + M x ^20+Mi 

(2) N 2 + M 2 ^ 2N + M 2 

(3) N 2 + N ^ 2N + N 

(4) i\TO + M 3 ^iV + 0-hM3 

(5) JV'O-hO^Oa + iV' 

(6) N 2 + O^NO + N 

(7) iVO ^ NO* 4* e~ (20) 

where Mi, M 2 , M3, are catalytic third bodies. The 
clean-air chemical model has 7 species (n = 7) and 
seven reactions (m = 7). In order to simulate the 
seeded air flow in the MHD channel, the potassium 
seeding reaction has been added to the above chem- 
istry model. This reaction is the ionization of atomic 
potassium (K) and is given by the following equation: 


(8) K+e~ ^K+ + e- + e~ (21) 

Using the law of mass action, the nondimensional 
mass production rate of species s is 


^ = M,fx, - u’ k<3 ) [K f , k (T) f[]jrr, 

fc=l r=X 

n t 


( 22 ) 


~K b , k (T) 
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where 7 r is the nondimensional mole-mass ratio of 
the reactants, M a is the molecular weight of species 
s, i/ ka and u k a are the stoichiometric coefficients and 
th is the number of reactants. Further details on 
the reaction rates and the thermodynamic and trans- 
port properties can be found in Ref. [26]. The electri- 
cal conductivity is determined from the species mole 
fractions, along with the temperature, density, and 
pressure of the gas, using the program of Park [21]. 


Numerical Method 

Solution of PNS Equations 

The governing PNS equations with MHD source 
terms have been incorporated into NASA’s upwind 
PNS (UPS) code [20]. These equations can be solved 
very efficiently using a single sweep of the flowfield 


for many applications. For cases where upstream 
(elliptic) effects awe important, the flowfield can be 
computed using multiple streamwise sweeps with ei- 
ther the IPNS [28], TIPNS [29], or FBEPNS [30] algo- 
rithms. This iterative process is continued until the 
solution is converged. 

For the iterative PNS (IPNS) method, the E vec- 
tor is split using the Vigneron parameter (cj) [31]. 
This parameter does not need to be changed for the 
present low magnetic Reynolds number formulation. 
In the previous high magnetic Reynolds number code 
[17] it was necessary to modify the Vigneron param- 
eter to account for MHD effects. After splitting, the 
E vector can be written as: 


E = E* + EP 


(23) 
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(24) 


The streamwise derivative of E is then differenced 
using a backward difference for E* and a forward dif- 
ference for the “elliptic” portion (E p ): 


(f|) = [(E;+1 ’ E:) + (E « ' Er " )1 

(25) 

where the subscript {i + 1) denotes the spatial index 
(in the f direction) where the solution is currently 
being computed. The vectors E* +1 and Ef +1 are then 
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linearized in the following manner: 


Solution of Species Continuity Equations 


E 


i+1 


E* + 



(U i+1 - Ui) 


E ? +1 = E ?+(^).(U. + 1 -U t ) 

The Jacobians can be represented by 


A* 
A f 


dE* 

dV 

dEP 

d\J 


(26) 


(27) 


After substituting the above linearizations into 
Eq. (25), the expression for the streamwise gradient 
of E becomes 


(f) M - 

+ (Ef+2 - Ef) ] (28) 

The final discretized form of the fluid flow equa- 
tions with MHD source terms is obtained by substi- 
tuting Eq. (28) into Eq. (16) along with the linearized 
expressions for the fluxes in the cross flow plane. The 
final expression becomes: 

[>*-«>+ 1 ( i ), ■ 

+ i(ro)r ,(Au ‘ ) ‘ +i=BHs(29) 

where 

(AU i ) k+1 = (U i+ i-U i ) fe+1 



and the superscript k+1 denotes the current itera- 
tion (i.e. sweep) level. In the preceding equation, the 
MHD source term, Smhd, is treated explicitly since it 
is evaluated using the velocity at station i (V*). For 
most cases, this will not degrade the accuracy of the 
solution since is small and the velocity changes 
slowly. If this is not the case, a predictor-corrector 
procedure can be implemented whereby a predicted 
velocity at station i+1 (Vjj^) is first obtained using 
Eq. (29). The solution at station i+1 is then recom- 
puted by evaluating Smhd with 


For chemical nonequilibrium, the species continuity 
equations, Eq. (19), must be solved in addition to 
the magnetogasdynamic equations. The equations 
have been integrated using the loosely-coupled ap- 
proach of Tannehill et ad. [32]. In this approach, 
the species continuity equations and magnetogasdy- 
namic equations are solved separately. The coupling 
between the two sets of equations is then obtained 
in an approximate manner. The species continuity 
equations are modeled using a second-order-accurate, 
upwind-based TVD scheme for the convective terms 
and second-order-accurate central differences for the 
diffusion terms. The assumption of zero net charge of 
the gas is used to eliminate the electron mass conser- 
vation equation. In addition, the species continuity 
equation for the nth species is eliminated by using 
the requirement that the mass fractions must sum 
to unity. The term representing the rate of produc- 
tion of species, uj 3 , is treated as a source term, and is 
lagged to the previous marching level. 

The coupling between the fluids and the chemistry 
is performed in an approximate manner. First, a fluid 
step is taken from marching station i to i+1 assuming 
frozen chemistry. Then the fluid density and velocity 
at i+1 are used in the solution of the species continu- 
ity equations to obtain species mass fractions at i+1. 
Finally, the species mass fractions, molecular weight 
of mixture, fluid density, and internal energy at i+1 
are used to obtain the new temperature, pressure, 
specific enthalpy, and frozen specific heats at the i+1 
marching station. 

The temperature is obtained by performing a 
Newton-Raphson iteration of the following form: 


Tk+ i_ T k g(T*)-ft 

am 


(30) 


where 

g( T) = J2c 3 h 9 (T) 

5=1 

5 '(t) = y> 3 c p , 4 (T) 

5=1 

and k is the index of iteration. The iterations are 
continued until 


| T k+1 - T k j < € 

where e is a small positive quantity. Once the temper- 
ature is determined, the pressure can be computed us- 
ing Dalton’s law of partial pressures. Further details 
of this procedure can be found in Refe. [32] and [33]. 
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Numerical Results 

The numerical calculation of the 3-D supersonic flow 
in the experimental MHD c h annel is now discussed. 
The flow in the nozzle section was computed using a 
combination of the OVERFLOW code [34] and the 
present PNS code (without MHD effects). For the 
OVERFLOW nozzle calculation, a highly stretched 
grid consisting of 150 x 80 x 80 grid points was used. 
The normal grid spacing at the wall was 1.0 x 10“ 5 
m. For the PNS calculation of the flow in the re- 
mainder of the nozzle and the rest of the MHD chan- 
nel, a highly stretched grid consisting of 90 points in 
both the y and z directions was used and the nor- 
mal grid spacing at the wall was 2.0 x 10~ 6 m. As 
a consequence of flow symmetry, only one-fourth of 
the channel cross section was computed in the 3-D 
calculations. 

The calculations were performed assuming turbu- 
lent flow throughout the MHD channel. The chan- 
nel wall temperature was assumed to be isothermal 
since quasi-steady flow conditions were maintained 
in the experiment for only about 1.2 milliseconds. A 
schematic of the powered portion of the MHD chan- 
nel along with the directions of the applied magnetic 
and electric fields is shown in Fig. 2. The values of the 
magnetic field ( B z ), and the electric field (E v ) were 
kept constant in the powered portion of the channel. 

Three different chemistry models were used in this 
study to simulate the flow in the MHD channel. 
These were: (1) perfect gas (7 = 1.25), (2) equi- 
librium air, and (3) nonequilibrium seeded-air chem- 
istry. For the perfect gas and equilibrium air calcu- 
lations, the electrical conductivity (cr e ) was assumed 
constant. For the nonequilibrium seeded-air calcu- 
lations, the electrical conductivity varied throughout 
the flowfield and was determined using the program 
of Park [21]. The seeding (as in the experiment) con- 
sisted of 1% (by mass) of potassium. Two test cases 
corresponding to Runs 15 and 16 of the NASA Ames 
experiments [2,35] were computed in this study and 
are now discussed. 


Test Case Is 

NASA Ames MHD Rim 15 

C Kap . = 320V) 

The dimensional 
are: 

flow parameters for this test case 

Vo = 

9.10 x 10 s N/m 2 

T 0 - 

5560 K 

T w - 

300 K 

<j e = 

130 mho/m (or variable) 

B z = 

0.0, 0.92 T 

Ey = 

0, 3955, 5000 V/m 


where the subscript o denotes total conditions at the 
nozzle entrance and w denotes wall conditions. 

This case was computed using several different elec- 
tric field strengths in order to properly simulate the 
experiment. In the experiment, the voltage applied 
to the electrodes was approximately 134 V for this 
case, however, due to the sheath voltage drop, the 
actual voltage applied to the flow is smaller than the 
electrode voltage. The voltage drop was measured for 
the central inviscid core flow, and was approximately 
67 V [2]. Since the boundary layer is computed in the 
numerical solution, the applied electric field must be 
approximately the voltage drop across the electrodes 
minus the sheath voltage drop. Unfortunately, it is 
not a trivial task to measure the sheath voltage drop. 
Therefore, several different electric field strengths 
were chosen in the numerical calculations so that 
the corresponding voltage drop across the electrodes 
would be between 67 V and 134 V. The voltage drop 
of 67 V corresponds to E v = 3955 V/m and a voltage 
drop of 84.7 V corresponds to E y = 5000 V/m. 

The computed streamwise variation of static pres- 
sure for the nonequilibrium seeded- air calculations 
is shown in Fig. 3 for the different electric field 
strengths. The pressure variation with no electric 
field or magnetic field is denoted by E v = 0. The re- 
sults for E v = 5000 V/m are in excellent agreement 
with the experiment. The numerical results show an 
increase in static pressure as the electric field strength 
is increased. The computed streamwise variation of 
static pressure for the different chemistry models is 
shown in Fig. 4 for E v = 5000 V/m. The nonequi- 
librium seeded- air model gives the closest agreement 
with the experimental pressures. 

For the nonequilibrium seeded-air computations, 
the electrical conductivity was not constant but var- 
ied throughout the flowfield. The average conduc- 
tivity (averaged over the channel cross section) at 
the center of the powered portion of the channel 
(electrode pair 10) was found to be 130 mho/m for 
E y = 5000 V/m. This value of conductivity is within 
the range determined in the experiments and is the 
same constant value that was used for the perfect gas 
and equilibrium air computations. 

The computed streamwise variation of averaged ve- 
locity for the nonequilibrium seeded-air calculations 
is shown in Fig. 5. The velocities are averaged across 
the channel cross section and normalized using the 
entrance velocity to be consistent with the experi- 
ment. In the experiment, the velocities were obtained 
by measuring the voltage generated by the flow at 
the last electrode pair (19) which is unpowered. This 
procedure inherently involves an averaging of the ve- 
locity profile. The numerical results indicate an in- 
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crease in the averaged velocity of about 27% with 
E y = 5000 V/m and this agrees exactly with the ex- 
perimental value of 27%. 

The computed streamwise variation of averaged ve- 
locity for the different chemistry models is shown in 
Fig. 6 for E y = 5000 V/m. Both the equilibrium air 
and the nonequilibrium seeded-air chemistry mod- 
els give similar results. The centerline variation of 
static temperature for the nonequilibrium seeded-air 
model is given in Fig. 7 for the different electric field 
strengths 


Test Case 2: NASA Ames MHD Run 16 

C Vcap. . = 380V) 

The dimensional flow parameters for this test case 
are: 

Po = 9.92 x 10 5 N/m 2 
T 0 = 5560 K 
T w = 300 K 

cr e = 130 mho/m (or variable) 

B z = 0.0, 0.92 T 

E y = 0, 4309, 5300, 6000 V/m 

This test case was also computed using several dif- 
ferent electric field strengths in order to properly 
simulate the experiment. The electric field strength 
of 4309 V/m corresponds to the voltage drop mea- 
sured in the central inviscid core flow. The com- 
puted streamwise variation of static pressure for the 
nonequilibrium seeded-air calculations is shown in 
Fig. 8 for the different electric field strengths. The ex- 
perimental pressure variation agrees with the numeri- 
cal result with an electric field strength of 5300 V/m. 
The computed streamwise variation of static pressure 
for the different chemistry models is shown in Fig. 9 
for E y = 6000 V/m. The different chemistry models 
produce similar results. 

The averaged electrical conductivity at the center 
of the powered portion of the channel (electrode pair 
10) was found to be 142 mho/m for E y = 5300 V/m. 
This value is at the upper end of the range deter- 
mined in the experiments and is higher than the con- 
stant value of 130 mho/m used in the perfect gas and 
equilibrium air computations. 

The computed streamwise variation of averaged ve- 
locity for the nonequilibrium seeded-air calculations 
is shown in Fig. 10. The numerical results indicate 
an increase in the averaged velocity of about 38% for 
E v = 6000 V/m and this agrees closely with the ex- 
perimental value of 39%. The computed streamwise 


variation of averaged velocity for the different chem- 
istry models is shown in Fig. 11 for E y = 6000 V/m. 
Once again, the equilibrium air and the nonequilib- 
rium seeded- air models give similar results. The cen- 
terline variation of static temperature for the equi- 
librium model is shown in Fig. 12 for the different 
electric field strengths. 

Concluding Remarks 

In this study, a new 3-D parabolized Navier-Stokes 
algorithm with nonequilibrium seeded-air capability 
has been developed to efficiently compute MHD flows 
in the low magnetic Reynolds number regime. The 
new algorithm has been used to compute the flow 
in the NASA- Ames experimental MHD channel for 
Runs 15 and 16. The numerical results are in good 
agreement with the experimental results. 
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Figure 1: Cross section of NASA Ames MHD Channel 



Figure 2: Schematic of powered portion of MHD channel 
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Figure 4: Streamwise variation of static pressure for Run 15 ( E y = 5000 V/m) 
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Figure 6: Streamwise variation of averaged velocity for Run 15 (E y = 5000 V/m) 
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Figure 8: Streamwise variation of static pressure for Run 16 (seeded-air model) 
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Figure 10: Streamwise variation of averaged velocity for Run 16 (seeded-air model) 
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Figure 12: Streamwise variation of static temperature for Run 16 (equilibrium air model) 
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